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A simple model of irreversible aggregation under differential sedimentation of particles in a fluid 
is presented. The structure of the aggregates produced by this process is found to feed back on 
the dynamics in such a way as to stabilise both the exponents controlling the growth rate, and 
the fractal dimension of the clusters produced at readily predictable values. The aggregation of ice 
crystals to form snowflakes is considered as a potential application of the model. 
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I. INTRODUCTION 



Simple models of cluster-cluster aggregation have been 
the focus of a great deal of interest, particularly over the 
last two decades. The structure of aggregates formed 
through a variety of dominating mechanisms (eg. diffu- 
sion limited reaction limited and ballistic motion 
^]) have been studied through theoretical, experimental, 
and computational work. 

Another aggregation mechanism which is relevant to 
several physical systems is that of differential sedimen- 
tation. Particles with a range of size and/or shape will 
almost inevitably sediment through a fluid at different 
speeds under the influence of gravity, leading to colli- 
sions. If there is some mechanism by which the particles 
stick on contact then aggregates will be formed. An ex- 
ample of this kind of phenomenon is the aggregation of 
ice crystals in Cirrus clouds. Small 'pristine' ice parti- 
cles are formed at the top of the cloud, and proceed to 
fall through it, colliding with one another and sticking to 
produce aggregates (snowflakes). 

The aim of this paper is to provide a simple model for 
growth by differential sedimentation which captures the 
essential physics of the system in the inertial flow regime, 
and to consider its application to snowflake formation. It 
is divided into five main parts - a description of the model 
and the assumptions underlying it; details of computer 
simulations and the results obtained from them; a theory 
section which offers an argument to account for the be- 
haviour observed in the simulations; an investigation of 
the model's applicability to snowflake formation; and a 
concluding discussion. The simulation results and their 
comparison to real cloud data have been presented briefly 
in a separate letter 0. 



II. MODEL 

We focus on the dilute limit, where the mean free 
path between cluster-cluster collisions is large compared 
to the nearest neighbour distance between clusters. In 
this regime we can limit our interest to individual binary 
collision events, ignoring spatial correlation. As further 
simplifying approximations, we assume that clusters have 
random orientations which do not signiflcantly change 
during a close encounter, that collsion trajectories are 
undeflected by hydrodynamic interaction, and that any 
cluster-cluster contacts result in a permanent and rigid 
junction. 

In order to sample the collisions between clusters, we 
flrst formulate a rate of close approach. For any two 
clusters i, j with nominal radii (see below) and rj re- 



spectively and fall speeds Vi 



, the frequency with 



which their centres pass closer than a distance (r^ -I- rj ) 
is proportional to the total area over which trajectories 
yielding a close approach event are possible, and the rela- 
tive speed of the pair. This is illustrated in flguren] The 
rate constant for approach closer than centre-to-centre 
separation (r^ -I- rj) is therefore given by: 



(1) 



In our computer simulations the nominal radii are cho- 
sen to fully enclose each cluster and the close approach 
rate calculated above is exploited to preselect candidate 
collision events. Collisions are accurately sampled by se- 
quentially choosing pairs of clusters with probability pro- 
portional to Tij , checking each pair for collision along one 
randomly sampled close approach trajectory, and corre- 
spondingly joining that cluster pair if they do indeed col- 
lide. In the theoretical arguments presented in section 
four, we make the simplifying assumption that all close 
approaches lead to collisions (or at least a fixed fraction 
of them do), using nominal radii based on fractal scaling 



from the cluster masses. 

The model is completed by an explicit form for the fall 
speeds entering equation We assume that the clus- 
ters are at most only partially penetrated by the fluid flow 
past them, so that cluster radius is the relevant length 
governing the drag force law. Then qualitatively and 
by dimensional argument we expect the same drag be- 
haviour as for a falling sphere, which may be written in 
the form: 



Fd = pVkf{Re) 



(2) 



where / is a function of the Reynolds number — rv /vk 
alone, p is the density of the surrounding fluid, and V}^ is 
the kinematic viscosity. Although details of the function 
f{Re) should be different from spheres, we still expect to 
have incrtial and Stokes regimes where / takes the forms: 



fiRe) 



i?g for incrtial flow 
Rp for viscous flow 



(3) 




fall speed v. 



Total area for 
all possible 
close approaches 



= TZ(r.+ r.) 

' J 



r. 

. J 




We consider below the general form f{Re) ^ Re^", with 
a as an adjustable parameter in order gain understand- 
ing spanning the two extreme regimes. Setting the drag 
force equal to the weight mg of the cluster, the terminal 
velocity is then given by 



FIG. 1: Illustration showing a possible scenario in which the 
centres of a pair of clusters falling at a relative speed \vi — Vj\ 
come within a distance {n + rj) of one another (a close ap- 
proach). The shaded circle illustrates the total area encom- 
passing all possible close approach trajectories = n{ri + rj)^ . 



r \pvl) 



(4) 



where a = ^ for inertial flow and a=l for viscous flow. 
A more complete discussion of the fall speed is given by 
Mitchell 0, but provided d/ > 2 so that cluster pro- 
jected area scales as the square of cluster radius, this 
reduces to a simple crossover between the above limits. 
The empirical crossover is very slow, spread over some 
three decades of Reynolds number, so fixed intermediate 
values of a can reasonably approximate behaviour over 
a significant range 0. In our simulations we took the 
radius determining the fall velocity to be proportional to 
the radius of gyration, and in our theoretical calculations 
we simply used the same nominal radii as for the collision 
cross sections above. 



III. COMPUTER SIMULATIONS AND 
RESULTS 

The primary particles used at the beginning of the sim- 
ulations were rods of zero thickness, half of which had a 
length (and mass) of unity, and half of which were twice 
as long and massive. Purely monodisperse initial condi- 



tions are not possible in this model, since \vi 



would 



be zero. Apart from this special case however, it is an- 
ticipated that the asymptotic behaviour of the system 
should be insensitive to the initial distribution, and in- 
deed the results described in this section appear to be 
preserved for a variety of starting conditions. 



In aggregation models it is typically the case (eg. Vic- 
sek and Family 0) that after the distribution has had 
time to 'forget' its initial conditions it will approach a 
universal shape. This is usually expressed by the 'dy- 
namical scaling' ansatz, which states that as m, s — s- oo: 



nm(t) = s{t) ^(j> 



m 



(5) 



where nm{t) is the number of clusters of mass m at 
time and the rescalcd distribution is a function of 
X = m/s{t) alone. The quantity s(t) is a characteris- 
tic cluster mass, and for non-gelling systems one expects 
that a suitable choice is given by the weight average clus- 



ter mass, s{t) = Ej'^I/Ej 



Using this choice our 



simulation data conform well to scaling, as shown in the 
left panel of figure |21 

The shape of the rescaled distribution was studied. A 
plot of (f>{x')dx' as a function of x is shown in the 
right panel of figure|21and shows an exponential decay for 
very large x, with a 'super-exponential' behaviour taking 
over as x approaches unity from above. This behaviour 
appears to be universal for all values of a in the range 
studied. 

For X <C 1 the qualitative form of (j){x) was found to fall 
into two distinct catagories depending on the value of a. 
For a > i the distribution appears to diverge as a power 
law: (l){x — > 0) shown in figure[5]for a — 0.55. 

The exponent r was found to be approximately constant 
at r ~ 1.6 ±0.1 over the range 5 < a < |. For a < 5 the 
distribution was found to be peaked, with a maximum at 
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0.8 



FIG. 2: Scaling of the cluster mass distribution. The left panel shows how the rescaled cluster size distribution (j) = s{t)^nm{t) 
converges to a universal function of rescaled cluster size x = m/s{t), where the data are overlayed for different values of the 
weight average cluster size, s{t) = 20, 50, 150, 400. The scales are logarithmic and a least squared fit (f>{x) ~ x~^'^ for x < 10~^ 
is shown by the dashed line. In the right hand panel J"" (l){x')dx' is shown on a semi-log plot, illustrating the exponential tail 
(dashed line is intended to guide the eye). Both simulations began with 250,000 rods, and used a = 0.55 in the sedimentation 
law. 



some small followed by a power law decay for 

< X < 1. 

Comparison with other aggregation models suggests 
that the clusters produced are likely to be fractal in their 
geometry, and in particular cluster mass and (average) 
radius should be in a power law relationship m ^ r'^f 
where d/ is the fractal dimension. A log plot of radius of 
gyration against mass for all the clusters produced over 
the course of the simulation is shown in figure |21 Also 
shown in this figure is the logarithmic derivative of the 
above plot, which shows the variation in the apparent 
fractal dimension of the clusters with size. From this 
plot, it seems that the fractal dimension approaches an 
asymptotic value as m ^ oo; in the case shown (a = 
0.55) we estimate this value as df ~ 2.2 ± 0.1. The value 
of this limiting fractal dimension was found to vary with 
a as shown in figure^ Note that our assumption df > 2, 
required to support the assumed fall speed relationship, 
is indeed satisfied for the physical range a> 1/2. 



IV. THEORY 

The most common theory used to describe cluster- 
cluster aggregation problems is that of von Smoluchowksi 
, which provides a set of mean- field rate equations for 
the evolution of the cluster mass distribution: 



dt 



1 °° 
= - ^ K,jn,{t)nj{t)-nk{t)^K^jnj{t). (6) 



clusters expressed (only) in terms of their masses i and 
J. Analytical solutions of Smoluchowski's equations have 
not been obtained except for a few special cases of Kij . 
However, Van Dongen and Ernst have shown that for 
non-gelling kernels (see below) the solutions approach the 
dynamical scaling form of equation ^ in the large-mass, 
large-time limit; substituting this into equations lO al- 
lows one to obtain some information about the asymp- 
totic behaviour of the rescaled cluster size distribution 

To apply this theory we need to compute the reaction 
rates Kij , which means averaging collision rate with re- 
spect to cluster geometry at fixed masses. This we es- 
timate by substituting averages from fractal scaling for 
the radii in equations ^ for the close approach rate and 
Q for the fall speeds, and assuming constant collision 
efficiency leading to: 



0' 



/df 



+ 3 



(7) 



Van Dongen and Ernst's analysis is sensitive to two expo- 
nents characterising the scaling of the coagulation kernel 
in the limit 1 <^ i ^ j, 



■LL -l 



which in our case yields: 

/I = min(0, a — dj^) 

v — max{a + dj^ ,2dJ^) 



(8) 

(9) 
(10) 



where rik (t) is the number of clusters of mass k at time 
t (per unit volume). The kernel Kij contains the physics 
of the problem, being a symmetric matrix, the elements 
of which govern the rate of aggregation between pairs of 



A third exponent combination X = ^ + u = a + dJ con- 
trols the growth of the average cluster mass through the 
differential equation s(t) = ws{t)^ , where w is a constant, 
and for the non-gelling case we require A < 1. 
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FIG. 3: Left hand panel shows a log plot of radius of gyration as a function of cluster mass for a — 0.55, averaged over four 
runs of 250,000 initial rods. Solid line indicates the theoretical prediction for the fractal dimension. The right hand panel shows 
the inferred fractal dimension as a function of cluster mass. Error bars are one standard deviation. Data points with a > 0.3 
have r 
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FIG. 4: Variation of the fractal dimension as a function of the 
parameter a. Circles are simulation data, solid line indicates 
theoretical prediction. 



Our identification of the exponent v is crucial to a 
mechanism by which the fractal dimension can control 
the dynamics. If the fractal dimension is low enough, 
then the exponent v will exceed unity. However, Van 
Dongen ^ has shown that the Smoluchowski equations 
predict the formation of an infinite cluster instantly in 
such a situation. In a finite system this clearly can- 
not occur, and it simply means that a few clusters will 
quickly become much larger than the others with their 
growth dominated by accretion of small ones. In this sce- 
nario the growth of the large clusters approaches that of 
ballistic particle-cluster aggregation, where it has been 
shown by Ball and Witten that the fractal dimen- 
sion of the clusters produced tends to df = 3. This in- 



creased fractal dimension reduces the value of v, forcing 
it back to a value of one if a < | . Through this feedback 
mechanism, a bound is placed on the fractal dimension 
df > max[2, (1 - a)~^] for a < |. 

The system could perhaps settle in a state where v <l. 
However, the growth in such a regime is much less biased 
towards collisions between clusters of disparate sizes, and 
the distribution is relatively monodisperse. This would 
tend to make collisions between clusters of a similar size 
likely, leading to much more open structures, with a lower 
fractal dimension, in turn acting as a feedback mecha- 
nism to increase the value of v. The authors suggest 
that, at least over some range of a, this effect will force 
the system towards the v = \ state. The discontinuity 
in the polydispersity of the system at v = 1 forces the 
system to organise itself such that it can remain at that 
point. This is similar to the argument put forward by 
Ball et al for reaction limited aggregation. 

If it is accepted that ;/ — > 1 then the fractal dimension 
of the clusters produced ought to be directly predictable 
from equation l(Tn|) : 

2 

d/ = max[2, (1 - a < -. (11) 

3 

A curve showing this theoretical behaviour is superim- 
posed on the simulation data in figure 01 and appears to 
show good agreement up to a ~ |. For a > | the theo- 
retical prediction is that df — 3 and v — a + ^ > 1, but 
because of its somewhat pathological nature we have not 
attempted to make simulations in this regime. It is how- 
ever clear from the extrapolation of our results in figure 
01 that this is likely to hold. 

Obtaining an exact form for the cluster size distribu- 
tion (t>{x) is a non-trivial exercise. However, following the 
methodology of Van Dongen and Ernst jg, we consider 



the small- a; behaviour of (/'(x) when df < (ie. < 0). 
In such a regime the small-x behaviour is dominated by 
collisions between clusters of disparate sizes; the gain 
term in the Smoluchowski equations may therefore be ne- 
glected, and one attempts to solve the integro-differential 
equation: w[x(j)' {x) + 2(f>{x)] = (f>{x) K{x,y)(f>{y)dy. 
For X <^ y, the kernel Q) may be approximated to 



K{x,y)^x^'y'' 



y , and one obtains: 



exp 



xf^p^ 



(12) 



where pi is the z"* moment of the rescaled distribution 
(f>{x), and the exponent r is given by r = 2 + pxw^^. 
It is clear that limj;^o[0(2;)] =0. As a; increases from 
zero, 4'{x) also increases, until reaching a maximum at 
Xm. = {wr/p^y/'-^ For Xm ^ a; ^ 1 the distribution 
has an approximately algebraic decay (f>{x) ~ x~'^ . This 
'bell-shaped' curve is consistent with the behaviour seen 
in the computer simulations when a < i. 

In the case df > a~^, it has been shown jQj that for 
all kernels with /i = 0, v < 1 the cluster size distribution 
diverges as a; with the form 



(13) 



where t = 2 — p\w ^ . This behaviour is consistent with 
the simulation for a > ^. The change in the qualitative 
shape of (f>{x <C 1) around a = then is further evidence 
to suggest that the system selects to sit at = 1. 

The shape of (j>{x 3> 1) has also been studied by Van 
Dongen and Ernst (ILj . They have shown that for non- 
gelling kernels, the tail of the distribution is expected to 
take the form 



(x) ^ a; 



(14) 



where 9 and S are constants. This would appear to be 
consistent with the behaviour observed in the simulations 
for all values of a, providing an exponentially dominated 
cut-off at large x. 



V. APPLICATION TO SNOWFLAKE 
FORMATION 

The principle motivation for the model presented in 
this paper was to attempt to understand some of the 
properties of Cirrus clouds. These are high altitude 
clouds with a base betwen 5,500 and 14,000 metres and 
they are usually composed solely of ice crystals [T^ . 
Amongst others, Heymsfield and Piatt have observed 
that the ice crystals in these clouds are predominantly 
composed of columns, bullets, bullet-rosettes and aggre- 
gates of these crystal types. It is these aggregates which 
we hope to model, since the dominant mechanism by 
which they grow is believed to be through differential 
sedimentation (eg. Field and Heymsfield 14]). We there- 
fore ignore the effects of diffusional growth, turbulence. 



a) EXPERIMENTAL IMAGES ^ #" ^ | 



b) SIMULATED AGGREGATES 

#4^ 



FIG. 5: Projected images of a) ice crystal aggregates ob- 
tained using a cloud particle imager (CPI, SPEC Inc., USA) 
during an aircraft flight through a Cirrus cloud at tempera- 
tures between — 44°C to — 47°C, and b) sample clusters from 
our simulations. 



mixing, and particle breakup, in order to concentrate on 
the effects of this mechanism alone. The Reynolds num- 
ber for aggregates of a few crystals is typically between 
~ 10 — 100 which ought to be modelled acceptably by our 
inertial flow approximation. Because we have not mod- 
elled the detailed hydrodynamics we may also be ignoring 
subtleties such as wake capture. 

All of the results below are presented for the purely 
inertial regime (assumed to be the most relevant to this 
problem) where a = ^. The initial particles were rods 
of zero thickness - however, the asymptotic behaviour is 
anticipated to be insensitive to the initial conditions, and 
indeed by running the simulation with 'bullet rosettes' for 
the initial particles (three rods, crossing one another at 
right angles, through a common centre), no change in the 
end results were found, only in the approach to scaling. 

Ice crystal aggregates have been studied through the 
use of cloud particle imagers during aircraft flights 
through ice clouds. Sample images from such a flight are 
shown in figure |31 alongside some of our simulation clus- 
ters. Using this experimental data, the geometry and size 
distribution of ice particles in these clouds has been stud- 
ied, allowing for quantitative comparison between theory 
and experiment. 

The fractal dimension of snowflakes in Cirrus may be 
inferred from the work of Heymsfield et al jl5i |. By 
measuring the effective density pe of bullet and bullet- 
rosette aggregates as a function of their maximum lin- 
ear dimension D, and fitting a power law to their data, 
they found the relationship pe ^ D~^-^^. This scaling 
implies that the aggregates have a fractal dimension of 
df = 2.04, which is consistent with the values predicted 
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FIG. 6; Aspect ratio of simulation clusters as a function of 
their maximum dimension. The curve seems to approach an 
asymptotic value of ~ 0.65, independent of the initial condi- 
tions used: here we show data for both rods and rosettes. 

by our model (simulation giving df = 2.05 ± 0.1 and the- 
ory giving d/ = 2). 

The aspect ratio of the clusters may also be calcu- 
lated. Random projections of simulation clusters were 
taken. The maximum dimension of the projection D was 
measured, as was the maximum dimension in the direc- 
tion perpendicular to that longest axis, Dy^. The ratio 
of these two spans were binned by maximum dimension, 
averaged, and plotted as a function of D as shown in fig- 
ure El The ratio quickly approaches an asymptotic value 
of approximately 0.65 ± 0.05. This compares well to the 
measurements of Korolev and Isaac [T^ , where the ratio 
seems to approach a value of ~ 0.6 — 0.7. 

Finally, the shape of the snowflake distribution of lin- 
ear size may also be compared with experiment. Field 
and Heymsfield presented particle size distributions 
of the maximum length D of ice particles in a Cirrus 
cloud. The data were obtained with an aircraft and repre- 
sent in-cloud averages of particle size distributions (num- 
ber per unit volume per particle size bin width) along 
15km flight tracks ranging from an altitude of 9500m 
(-50°C) to 6600m (-28°C). To compare this data to 
the distributions obtained from simulation, we first nor- 
malise the data, and then make use of the dynamical 
scaling form jSJ, to collapse the distributions onto a sin- 
gle curve. Details of this are given in the appendix to 
this paper. The resulting histograms are shown in figure 
Hand appear to show quite good agreement, given the 
level of approximation present in our model. 



VI. DISCUSSION AND CONCLUSIONS 

A simple mean-field model of aggregation by differen- 
tial sedimentation of particles in an inertial flow regime 



has been constructed, simulated by computer, and anal- 
ysed theoretically in terms of the Smoluchowski equa- 
tions. It has been shown that there is strong numerical 
evidence, in addition to a theoretical argument, to back 
up the idea that the polydispersity of the distribution and 
the fractal dimension feed back on one another in such a 
way as to stabilise the system at v = 1. Above this value, 
the dominance of collisions between clusters of very dif- 
ferent sizes is so great as to push df towards a value of 
three. This in turn pulls the exponent v back down to 
unity. For v < 1 the system is quite monodisperse, re- 
sulting in relatively many collisions between clusters of 
similar sizes, and the fractal dimension is reduced, forcing 
V back up. The discontinuity in the shape of the distribu- 
tion around v = 1 is thought to provide the mechanism 
by which the system can stabilise at that point. 

If it is accepted that — > 1 , then the fractal dimension 
of the clusters produced may be predicted, and figure 
^ shows that this prediction agrees well with simulation 
results for < a < |. The sudden change in the be- 
haviour of df{a) and in the small- cc form of the cluster 
size distribution around a = ^ is strong evidence for the 
self-organisation proposed between df and v. 

For a > I the system is forced into a regime where 
v > 1, which has been regarded as unphysical because 
the Smoluchowski equation ^ predicts infinite clusters 
in zero time 0]. In the light of our results this regime 
merits further study beyond the Smoluchowski equation 
approximation (17i |. The value a = 1 is given by viscous 
flow, but here our form for does not include all of 
the relevant physics: in particular, small clusters may be 
caught in the fluid flow, and swept around larger clusters 
rather than hitting them, reducing the dominance of big- 
little collisions. This has been discussed in more detail 
for the particle-cluster aggregation case by Warren et al 

m 

The application of the model to the formation of ice 
crystal aggregates in Cirrus clouds has been considered: 
the fractal dimension, aspect ratio, and shape of the clus- 
ter size distribution seen in the model were all found to be 
consistent with experimental data. This is a promising 
indication that the ideas presented in this paper may be 
an acceptable model for the essential physics of snowflake 
aggregation in Cirrus. 
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SIMULATION EXPERIMENT 
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FIG. 7: The left hand panel shows the distribution of clusters by linear size at various stages of the simulation, rescaled in such 
a way as to collapse the data (see appendix). Initial conditions were 250,000 rods and the parameter a was set to a value of i. 
The right hand panel is a test of the same scaling using the experimental data presented by Field and Heymsfield [l^ . 



APPENDIX: SCALING OF THE CLUSTER 
RADIUS DISTRIBUTION 

Experiments have reported the distribution of ice ag- 
gregates by linear span rather than by mass, and we 
present here how that distribution ^ should naturally 
be rescaled. This tests the dynamical scaling ansatz 
which, for the mass distribution, gave ^ = rim = 
s~°'(j){m/ s), where a = 2 in mass-conserving systems. We 
anticipate fractal scaling so that m ^ D'^f and hence: 

From this expression we may calculate the moments of 
the distribution M{b) = J ^DMD in terms of the av- 
erage cluster mass s{t): 

M(6) - s-''+i+''/'*W x''/'^f(f>{x)dx (A.2) 

Jl/s 

where x = m/s. At small sizes we expect (t>{x) x^'^ . If 
b > rf/(r — 1) therefore, the integral converges as s — > oo. 



and M(6) ^ g-a+i+b/df ^ From our simulations, we have 
measured t~1.6, c?/~2, and so the lowest integer mo- 
ment which scales in this way is the second. We therefore 
choose this to normalise our data: 



[Mi2)]-'— D^f-'s-'-^^^f(b(-) (A.3) 
dD \ s / 



which, defining the average cluster diameter D* = 
M(3)/M(2) - s^Z-^f yields: 



[M{2)r^^{D*)-'^(^-§,y (A.4) 



where = Hence, if we as- 

sume that df approaches a constant value, plots of 
{[M{2)]-^.^.{D*)^} against (D/D*) should all lie on 
a single curve. 
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